Improved timing formula for the PSR B1257+12 planetary system 
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ABSTRACT 



We present a new analysis of the dynamics of the planetary system around the 



pulsar B1257+12. A semi-analytical theory of perturbation between terrestrial-mass 
planets B and C is developed and applied to improve multi-orbit timing formula for 
this object. We use numerical simulations of the pulse arrival times for PSR B1257+12 



to demonstrate that our new timing model can serve as a toll to determine the masses 
of the two planets. 



O 

Subject headings: celestial mechanics — planetary systems — pulsars: individual (PSR 
B1257+12) 



1. Introduction 

The first extra-solar planetary system has been discovered around a millisecond pulsar B1257+12 
(Wolszczan &: Frail 1992). The system consists of three planets named A,B and C with planets B 
and C having the orbits close to a 3:2 commensurability. This circumstance allows us to analyze 
the dynamics of the system beyond the classical Keplerian approximation. Namely, in such config- 
uration, the gravitational interactions of planets B and C give rise to observable time variations of 
B and C orbital elements. It was thoroughly discussed by Rasio et al. (1992) and Malhotra (1993a) 
(see also Malhotra et al. 1992; Malhotra 1993b; Rasio et al. 1993; Peale 1993). Subsequently, 
these studies were used to confirm the existence of the PSR B1257+12 planetary system through 
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detecting, in the timing observations of the pulsar, the presence of those non-keplerian variations 
(Wolszczan 1994) 

However, the application of non-keplerian dynamics goes further than the confirmation of the 
discovery. It can be used to derive some interesting information about the system which is not 
otherwise accessible. The aim of this paper is to apply the theory of perturbed planetary motion 
and derive an improved model for the timing observations of PSR B1257+12. Such model, as we 
show on simulations, should lead to determination of the masses of planets B and C, as well as 
the inclinations of their orbits. To this end, in section 2 we analyze the equations of motion in 
the barycentric and Jacobi coordinate systems, which we use in the paper. In section 3 we show 
how these two slightly different representations are related to the commonly used timing model for 
pulsars with companions. In section 4 we demonstrate how to express such timing formula in terms 
of the osculating orbital elements. In section 5 we show how to obtain the osculating elements of 
B and C. In section 6 we present an improved timing model describing the motion of this system 
and finally, in section 7, we perform numerical tests which show how it can be used in practice. 



2. Equations of Motion 

Let us consider a system consisting of a neutron star Pq with the mass mo, and N planets Pj 
with masses rrii. In an arbitrary inertial reference frame equations of motion of this system have 
the form 



N 

miirij 

3=0 %3 
3+i 



m l R, = -G^^(R i -R J ), (1) 



where 

Rij = URi-Rj-ll, i,j = 0,...,N, 
and G is the gravitational constant. The Hamiltonian function for system (1) has the form 

i=0 1 0<i<j<N lJ 

and equations (1) can be written as Hamilton's equations 

d n dH d „ dH 

Analytical perturbation theory for a planetary system is usually formulated in the so-called Jacobi 
coordinates r^, i = 0, . . . , N which are defined in the following way 

j k-l 

r k = R fc m iRi, for k = l,...,N, (4) 
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and 

N k 



ro 

UN 



— y^mjRj, fi k = y^mi. (5) 

//. AT ' * ' * 

i=Q i=0 

In other words, rj is the radius vector from the center of mass of bodies Pq, . . . , to body Pi, 
and ro is the center of mass of the system. The above formulae define a one-to-one relationship 
between Cartesian inertial coordinates and Jacobi coordinates. The inverse relationship has the 
form 

N 

K k = r + ^v k - — r - k = 0,...,N, (6) 
i=T + i * 

where we assume /it_i = 0. In terms of canonical Jacobi coordinates, Hamiltonian (2) reads 

Mat ^IH-irrH 1 ^ n (7) 
+ #i (ri, . . . ,r N ) , 

where H\ = i?i(ri, . . . , rjy) is its perturbative part. If we assume that masses of planets are small 
and all are of the same order e, then H 1 is of order e 2 , i.e., H\ = Hi + 0(e 3 ), where 



Hi = —G mimj 

l<i<j<N 



1 r ^ • rj - 



r, 



(8) 



From the form of Hamiltonian (7) it follows that po is a first integral and that the equations for 
variables {ri, ...,rjv,pi,...,pjv} do not depend on its particular value. Thus, we can assume that 
Po = implying that ro is constant and can be set to ro = 0. This is equivalent to the assumption 
that our inertial frame is a certain barycentric reference frame. The dynamics of the system is 
governed by Hamilton's equations with the Hamiltonian of the form 

H = H + H 1 + ---, (9) 

where 

N N 

tf = £^P?-G£«, (io) 

i=\ i=\ 

is the unperturbed part of the Hamiltonian describing a system of N independent planets. It 
follows from the form of Eq. (10) that each planet moves in a Keplerian orbit in the same way as a 
body with the mass m^-i/^ around a fixed gravitational center m^mi. Each of these Keplerian 
motions can be parameterized by Keplerian elements {T p , a, e, i, u, £1} — the time of pericenter, semi- 
major axis, eccentricity, inclination, argument of pericenter and the longitude of ascending node, 
respectively. 
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3. Timing model and coordinate systems 

Timing observations of pulsars represent measurements of the times of arrival of pulsars pulses 
(TOAs). An extraordinary precision of timing measurements allows a detection of very low-level 
effects in timing residuals (see for review Lyne & Graham-Smith 1998). In the case of a binary 
pulsar the observed TOAs exhibit periodic variations resulting from the motion of the pulsar around 
the center of mass. Such variations are modeled with the formula 



At = x 

where 



(cos E — e) sin oj + \A — e 2 sin E cos oj , (11) 

x = asmi/c, E — e sin E = n (t — T p ) , n = ~pi 

and At is the additional delay/advance in TOAs, E is eccentric anomaly, a, e, u, smi, P, T p are 
Keplerian elements of the orbit and c is the speed of light. When a pulsar has N planets the TOA 
variations become 

N 

At = \ Xj (cos Ej — ej) sin u>j + 

U L (12) 

+ yj\ — e 2 j sin Ej cos ojj) 

Note that in practice the parameters of the model which we determine, by means of the least-squares 
fit to the data, are Xj,ej,ujj, Pj,T p j, i.e., the projection of semi-major axis of the pulsar's orbit, 
eccentricity, argument of pericenter, orbital period and time of pericenter. In order to precisely 
understand and interpret these parameters we describe the pulsar's motion in the barycentric 
reference frame with the z-axis of the system directed toward the barycenter of the solar system 
and xy plane of the reference frame in the plane of the sky. This way, we have 

1 1 N 

At = — R -Z, R = VmA, (13) 

c to t-f 

where Z is the unit vector along the 2-th axis of the pulsar system barycentric frame and Rj are 
barycentric positions of planets. Assuming that 

m 'i ■ ■ G m 2 3 /-. a\ 

dj smij = XjC, — ; = n.-a,-, (14) 

to (1 + mj/moY J J 

with planets ' orbital parameters aj , rij , ej , ujj , T p j , we can most naturally interpret the motion of 
the pulsar as a superposition of the elliptic motions of its planets around the barycenter of the 
system. However, as it was mentioned in the previous section, the analytical perturbation theory 
is usually formulated in Jacobi coordinates in which the TOA variations become 

i N 
At = — R • Z, R = - 2^ KjTj, Kj = (15) 

C 3=1 ^ 
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where Vj are positions of planets. Furthermore, we have the following relations 

. . GmoHj Gm 2 3 . . 

K-jdi smij = XjC, = = n^CLA. lb) 

J J J J 1 - Kj J J 

Thus from the timing formula for TOA variations of a pulsar with planets it is possible to obtain 
two somewhat different descriptions of the pulsar motion. Although, they both represent a sum 
of a certain number of elliptic motions, the interpretation of some of their parameters is slightly- 
different. Throughout the rest of this paper, we will use Jacobi coordinates as they are more 
convenient in the formulation of the theory of perturbed motion. 



4. Osculating orbital elements 

The non-keplerian motion of the PSR B1257+12 system can be described by means of the 
osculating ellipses (i.e. by means of ellipses which parameters change with time). The time evolution 
of orbital elements can be determined by solving the classical Lagrange's perturbation equations 
(see, for example Brouwer & Clemence 1961) with the perturbation Hamiltonian given by equation 
(8). Such approach leads to the solution in the form 

x = x° + Ax(t-to), (17) 

where x stands for a specific orbital element, x° its initial value and Ax for its time dependent part 
of small magnitude Ax(i — io)/ x ° ^ 1- in the case of the PSR B1257+12 planetary systems the 
most significant part of the perturbations comes from planets B and C. Therefore, we can assume 
that the orbital elements of planet A are approximately constant while the elements of planets B 
and C change with time. Thus using the formulae for At from the previous section (Eq. (12)) and 
assuming the time evolution of orbital elements in the form (17), the additional TOA variations 
5tj,j = {B, C} due to the interactions between planets B and C can be expressed as follows 



cStj 



Kj sin i® 



where 



+ ^A^a°sin(2A°) + 

+ Aa,j sin X°j + A\ja°j cos A°, 



hj = ej sin uij , kj = ej cos uij , 



+ 

(18) 



(19) 

\j =nj(t-T pj )+ujj, 

and equation (18) is given to first order in Aaj, AXj, Ahj, Akj and the lowest order in ej (or hj 
and kj, since the eccentricities of planets B and C are very small). The above equations can be 
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obtained from the following well-known expansions 

co Sj B = -| + 2^ -J k _ 1 (ke)cos(kM), 

2 1 keZ ° (20) 
s'mE = - ^ -j^Jk~i{ke) sin(feM), M = A - to 
kez 

where J k (x) is a Bessel function of the first kind of order n and argument x, Zq denotes the set of 
all positive and negative integers excluding zero; and an approximate relation for Bessel functions 
of small arguments 

k 

Mx) ~ x « 1 (21) 

Now, our next step is to find the explicit form of the functions 

aj (t) = a° + Aa.it), X 3 (t) = A° + AA,(t), 

(22) 

hj (t) = h] + Ahj(t), kj (t) = fcj + Akj (t) . 



5. Semi-analytical perturbation theory 

Mutual gravitational interactions between planets cause periodic and secular changes of their 
orbital elements. In the case of PSR B1257+12 periodic variations can be related to conjunctions 
('close encounters') of planets B and C with the frequency n ce = rig — nc and to the 3:2 near- 
commensurability with the frequency n r = 3nc — 2ng. The dynamics of this system was studied 
by Rasio et al. (1992) and Malhotra (1993a) (see also Malhotra 1993b; Rasio et al. 1993; Peale 
1993). Malhotra (1993a) solved the Lagrange's perturbation equations for the orbital elements 
assuming non-resonant system with coplanar orbits. This first order perturbation theory is valid 
for (sini) -1 > 10 (orbital inclinations larger than about 6 degrees) which corresponds to mass 
ratios rrij/mo >6x 10~ 5 (and non-resonant case). When planets B and C are locked in the exact 
resonance it is necessary to develop a different theory of motion (Malhotra et al. 1992; Rasio et al. 
1993). 

It turns out that the first-order perturbation theory for this system developed by Malhotra 
(1993a) is not accurate enough for the purpose of determination of the planets' masses. Therefore, in 
this paper we present a new approach to this problem which precisely addresses the issue. Namely, 
first we assume that the orbits of B and C are not coplanar however their relative inclination I 
is small. The geometry of the system can be represented as in Fig. 1 and the perturbative part 
of Hamiltonian expanded to the first order in / and the product of the mass of planets B and C, 
?nim2, has the following form 



Gmim 2 
Hi = 

T2 



1 — 2— cos?/> + ( — ) I ——costp 
r2 \r2j / r 2 



(23) 
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where 

1p = fl + Ul - h ~ ^2 - T, T = T 1 -T 2 , 

a 3 (l - e 2 ) 

ri = Tf f ' J = 1.2, 

1 + ej cos /j 

and fj is the true anomaly. Note that in the above expansion, H\ does not depend on /. The first 
term in which / appears is of the second order; precisely it is sin 2 (//2). Thus as long as sin 2 (//2) 
is negligible, the above expansion is valid. 

Next, from the classical form of Lagrange's perturbation equations (see, for example Brouwer 
& Clemence 1961) we can obtain the following set of the first order differential equations for the 
elements a,j , ej , Uj , Xj 

-2 dHi 

3 i)n i ' ( 24 ) 



~ 3 dujj j 



A,- = rij + 



dHi 



l-e] 



dHi 



m a 3 <'>"., m a] (l + y/l - efj 



de-j 



where rij is the mean motion and <jj is related to time of pericenter through <jj = —nfT V j so as the 



Kepler equation reads 



ej sin Ej 



nd + (Tj = Xj 



(25) 



and the true anomaly /,• is related to the eccentric anomaly Ej through 



tan 



fj 



'l + ej 
1 — 



tan • 



(26) 



And it should be noted that the elements Qj , ij remain approximately constant in the case of small 
relative inclinations. 

In principle, such set of equations could be solved for aj,ej,Uj, Xj. In practice however, it is 
extremely complicated. On the other hand, in fact we do not need analytical formulae such as 
those presented in Malhotra (1993a). Thus the most suitable approach is to solve this problem 
numerically. In order to do so we just need the explicit form of the right-hand-side functions of 
equations (24) which can be easily obtained using the perturbing Hamiltonian H±. Subsequently, 
the equations can be solved numerically for a,j , ej , Uj , Xj . As we show in section 7, this approach 
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gives very accurate results, much more accurate than any reasonable analytical treatment. From 
the form of equations (24) it follows that the change of orbital parameters Ax (strictly speaking here 
Ax stands for AA, Ah, Ak and Aa/ao) is proportional to mi/mo and m2/mo. Precisely, corrections 
Axi are proportional to m-ij 'mo and Ax2 are proportional to m\ /tuq. Let us finally note that in 
this model there is one additional parameter r which in general is not known a priori. However, it 
can be determined through a least-squares analysis of the data. 



6. Timing formula 

We have now all the components necessary to derive a useful timing formula which will describe 
the motion of the PSR B1257+12 system. First of all, let us note that due to the small mass of 
planet A we will have: 

m A 

k a = Kl = , 

m 

(27) 

KB = «2 = ; ; ~ : , 

mo + "\4 + m B rn + m B 

mc mc 



KC = K 3 



mo + m.4 + m B + m c m + mg + mc 



The problem of finding masses and inclinations of planets B and C can be now formulated as 
a least-squares problem in which we fit to the data the following function 

At = At K (V) + 5t B (V, 7H , 7C ) + 8tc(V, 7s, 7 c), (28) 

where 

* = {x%, n° B , e%, u° B , T p ° e , x° c ,n° c , e° c , lo° c , T p ° c } , 

mg mc 

IB = , 7C = , 

m m 



-9- 



Atxi^) describes the Keplerian part of the motion given by equation (12) and 



StB = x° B 



+ 



+ 



^Ak B sm(2X B ) + 

Aa B . . n a \ ,n 
H — sm Ag + AAg cos X B 



Stc = x° c 



- Ah c ( ^ + J cos 



+ 



^ + Icos(2A°)) + 
^A/cc sin (2A C ) + 

+ sin A c + AA C cos A° 



(30) 



describe its non-keplerian part. In this model the initial values of osculating orbital elements replace 
the Keplerian elements as parameters of the model and we additionally have the parameters, ^ B 
and 7c, related to the masses of B and C. We also need the derivatives of At with respect to the 
parameters 7b, 7c}- With sufficient accuracy the derivatives with respect to ^ can be computed 
from the Keplerian part Atx of the model only and the derivatives with respect to 7^3 and 7c can 
be easily obtained as 5t B and Stc ar e proportional to them 



OAt 

dj B 



Stc 

IB ' 



dAt 

die 



Sts 
7c ' 



(31) 



Now, determination of masses of B and C and inclinations % and ic can be carried out in the 
following way. First, we assume that the mass of the pulsar mo is canonical (i.e. 1.4M Q ) and derive 
mg, mc directly from j B and 7c. Subsequently inclinations can be computed from the following 
equations 



KjCpj sin ij, 



Girio 



1 



o 2 o 3 



(32) 



where j = {B, C}. 



7. Numerical tests and conclusions 

We start the tests with the comparison of our approach for computing osculating elements with 
that of Malhotra (1993a). It turns out that the most significant component of the non-keplerian 
motion comes from the changes of mean longitude Xj. Therefore in Fig. 2 we present AXj for the 
case of coplanar, edge-on orbits. As we can see, the approach developed in this paper is in fact 
as good as the integration of full equations of motion. For other elements we obtain very similar 
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results. It should be mentioned that because the model of Malhotra (1993a) is not accurate enough 
in predicting the secular change of AXj and because AXj are proportional to 70 and 7c, using such 
model would result in a significant error in determination of the planets' masses. Next, we compared 
the TOA residuals due to the non-keplerian part of the motion calculated from our model with 
those obtained by means of numerical integration of full equations of motion. In Fig. 3, as one 
can see, for small relative inclinations (Fig. 3 (a) and (b)) our model is very accurate. For larger / 
(Fig. 3 (c) and (d)) one can see small differences but this is consistent with the assumptions made 
to obtain the model. 

Finally, we performed the tests which show how the derived timing model can be used to obtain 
the masses of planets. To this end we used the TEMPO software package (see the Internet location 
http : //pulsar .princeton. edu/tempo) modified to include our model of the non-keplerian motion 
of PSR B1257+12 . We simulated two different sets of artificial TOAs assuming the parameters of 
PSR B1257+12 (Wolszczan 1994). The first set of TOAs sampled every day covered a period of 
10 years. We also added Gaussian noise with a = 0.1/xs. The second set was prepared to resemble 
the real timing observations of PSR B1257+12. Under such conditions we simulated TOAs for the 
cases described in Fig. 3. Subsequently, we applied the modified TEMPO to obtain the masses of 
planets B and C. The results are presented in Table 1. 

The tests performed on Set 1 were used to establish a limit on applicability of our model. As 
one can see, for the relative inclinations / of about 10 degrees the model becomes inaccurate and 
the relative error of determination of masses is at the level of 20%. Because the relative inclination 
weakens the interactions of planets, using the model which assumes a small /, results, in such 
situations, in the determined masses that are smaller than the real ones. From the tests performed 
on Set 2 we also learn that this error is bigger than uncertainty originating in TOA measurement 
errors. On the other hand, when / is relatively small we obtain a very accurate determination of 
masses thus proving that our model can be successfully applied as long as the assumptions made 
to derive it are satisfied. 

To sum up, our model can be applied in all the cases when the 'observational' uncertainties 
of the planets' masses are bigger than the uncertainties resulting from the assumptions made to 
derive the timing model. This, in general, depends on masses of planets and the relative inclination 
of the orbits I as well as the characteristic of the observations but one can estimate that in the 
case of the PSR B1257+12 data / should be smaller than 10 degrees. We should also mention 
that for non-coplanar orbits the angle r will be in general different from zero therefore we have to 
find it in order to get the proper values of the planets' masses. It can be done by computing the 
least-squares value of x 2 for a range of r and then choosing the r that corresponds to the minimum 
of x 2 . And eventually, one has to remember that the PSR B1257+12 planetary system could be in 
such configuration that a more significant alteration must be done to our model. First of all, the 
relative inclination of the orbits could be large. In such case one would have to use the full form of 
the perturbative Hamiltonian H\ which would lead to a much more complicated model with four 
additional parameters Qq, Jlc, ig, ic (instead of one r). Secondly, the presence of a massive distant 
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planet, if significant for the theory of motion, would have to be taken into account. 
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Fig. 1. — Geometry of the system. The angles u±, L02 are the arguments of pericenter, Q±, 0,2 
longitudes of ascending node and ii, 12 inclinations of the orbit of the planets B and C; t\ and T2 
are the angles n\On\2 and ri20n\2 respectively. The angles I,ti,T2 can be found by solving the 
spherical triangle n\n\2'n2- 

Fig. 2. — Time evolution of the mean longitude AA(i) = X(t) — X°(t) for the planets B and C in the 
case of coplanar, edge-on orbits. The solid line indicates the solution by means of the numerical 
integration of full equations of motion, the open circles indicate the solution obtained with the 
approach presented in this paper and the dash-dotted line with the model by Malhotra (1993a). 

Fig. 3. — TOA residuals due to the non-keplerian part of motion of the PSR B1257+12 in four dif- 
ferent configurations. The solution obtained by means of the numerical integration of full equations 
of motion is indicated with the dots and the one obtained using our model with the solid line. 
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Table 1. Assumed and derived parameters of simulations a 



Assumed Set 1 Set 2 





T 


I 


m B 


mc 


m B 


m c 


m B 


mc 








[M e ] 


[M e ] 


[M ffi ] 


[M ffi ] 


[Me] 


[Me] 


(a) 


0?0 


0?0 


3.41 


2.83 


3.42(1) 


2.82(1) 


3.53(51) 


2.63(51) 


(b) 


0?0 


2?0 


4.99 


4.32 


4.98(1) 


4.26(1) 


5.07(51) 


4.06(48) 


(c) 


0?0 


10?0 


4.82 


4.94 


4.06(3) 


4.11(3) 


4.06(51) 


3.79(48) 


(d) 


9? 7 


10?3 


9.96 


16.33 


8.42(12) 


13.39(12) 


8.25(45) 


13.16(42) 



Figures in parentheses are 3a uncertainties in the last digits quoted. 
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Fig. 2.- 
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